# ============================================
# MODELO DEFINITIVO
# OMP (k fijo=4) + CV temporal (ventana fija + gap) + estabilidad
# Escala SOLO columnas float; dummies 0/1 (int32) pasan sin escalar
# X_train, X_test, y_train, y_test ya creados
# ============================================
import numpy as np
import pandas as pd
from sklearn.model_selection import GridSearchCV, TimeSeriesSplit
from sklearn.preprocessing import StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import OrthogonalMatchingPursuit
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
from sklearn.base import clone
from collections import Counter
# ---------- utilidades ----------
def rmse(y_true, y_pred):
return float(np.sqrt(mean_squared_error(y_true, y_pred)))
def adjusted_r2(r2, n, p):
if n - p - 1 <= 0:
return np.nan
return float(1.0 - (1.0 - r2) * (n - 1) / (n - p - 1))
def mape(y_true, y_pred):
y_true = np.asarray(y_true, dtype=float).ravel()
y_pred = np.asarray(y_pred, dtype=float).ravel()
denom = np.where(y_true == 0, np.nan, np.abs(y_true))
return float(np.nanmean(np.abs(y_true - y_pred) / denom) * 100.0)
def strip_remainder_prefix(names):
return [n.replace("remainder__", "") for n in names]
# ---------- preprocesado ----------
# Estandariza SOLO columnas float; resto (incl. dummies int32) passthrough
float_cols = list(X_train.select_dtypes(include=['float32','float64']).columns)
preprocess = ColumnTransformer(
transformers=[("num", StandardScaler(with_mean=True, with_std=True), float_cols)],
remainder="passthrough",
verbose_feature_names_out=False
)
# Pipeline base
pipeline = Pipeline(steps=[
("prep", preprocess),
("omp", OrthogonalMatchingPursuit(fit_intercept=True))
])
# ---------- CV temporal con misma ventana/gap ----------
VENTANA = 72
GAP = 12
try:
tscv = TimeSeriesSplit(n_splits=5, max_train_size=VENTANA, gap=GAP)
except TypeError:
tscv = TimeSeriesSplit(n_splits=5, max_train_size=VENTANA)
# ---------- Diagnóstico de k en {3,4,5,6} (sin 1-SE) ----------
param_grid = {"omp__n_nonzero_coefs": [3, 4, 5, 6]}
grid = GridSearchCV(
estimator=pipeline,
param_grid=param_grid,
scoring="neg_mean_absolute_error",
cv=tscv,
n_jobs=-1,
refit=True,
verbose=0,
return_train_score=True
)
grid.fit(X_train, y_train)
cv_results = pd.DataFrame(grid.cv_results_)
cv_results["k"] = cv_results["param_omp__n_nonzero_coefs"].astype(int)
cv_results["MAE_CV"] = -cv_results["mean_test_score"].astype(float)
print("Diagnóstico CV (misma ventana/gap):")
print(cv_results[["k","MAE_CV"]].sort_values("k").to_string(index=False))
# ---------- Modelo final para evaluar en test con k=4 ----------
K_FINAL = 4
model_k = clone(pipeline).set_params(omp__n_nonzero_coefs=K_FINAL)
model_k.fit(X_train, y_train)
# Predicciones y métricas
y_pred_train = model_k.predict(X_train)
y_pred_test = model_k.predict(X_test)
mae_test = float(mean_absolute_error(y_test, y_pred_test))
rmse_test = rmse(y_test, y_pred_test)
mape_test = mape(y_test, y_pred_test)
r2_tr = float(r2_score(y_train, y_pred_train))
r2_te = float(r2_score(y_test, y_pred_test))
# p = nº de coeficientes ≠ 0 (en el espacio transformado; válido para OMP)
coefs_std = model_k.named_steps["omp"].coef_
p = int(np.sum(np.abs(coefs_std) > 0))
r2_adj = adjusted_r2(r2_te, n=len(y_test), p=p)
metricas = pd.DataFrame({
"métrica": ["MAE test", "RMSE test", "MAPE test (%)", "R² test", "R² ajustado test", "R² train"],
"valor": [mae_test, rmse_test, mape_test, r2_te, r2_adj, r2_tr]
})
print(f"\nMétricas (modelo k={K_FINAL}):")
print(metricas.to_string(index=False))
# ---------- Coeficientes a escala ORIGINAL + % de importancia ----------
prep = model_k.named_steps["prep"]
feat_names_raw = list(prep.get_feature_names_out())
feat_names = strip_remainder_prefix(feat_names_raw)
mean_vec = np.zeros(len(feat_names_raw), dtype=float)
scale_vec = np.ones(len(feat_names_raw), dtype=float)
if len(float_cols) > 0:
scaler = prep.named_transformers_["num"]
raw_index = {name: i for i, name in enumerate(feat_names_raw)}
cont_idx = [raw_index[c] for c in float_cols]
mean_vec[cont_idx] = getattr(scaler, "mean_", np.zeros(len(cont_idx)))
scale_vec[cont_idx] = getattr(scaler, "scale_", np.ones(len(cont_idx)))
coefs_orig = coefs_std / scale_vec
b0_std = float(model_k.named_steps["omp"].intercept_)
intercepto_original = b0_std - float(np.dot(mean_vec / scale_vec, coefs_std))
mask = np.abs(coefs_std) > 0
sel_features = np.array(feat_names)[mask]
sel_coef_std = coefs_std[mask]
sel_coef_ori = coefs_orig[mask]
abs_sum = float(np.sum(np.abs(sel_coef_ori)))
porc = (np.abs(sel_coef_ori) / abs_sum * 100.0) if abs_sum > 0 else np.zeros_like(sel_coef_ori)
selected_features = pd.DataFrame({
"feature": sel_features,
"coef_estandarizado": sel_coef_std,
"coef_original": sel_coef_ori,
"porc_importancia": porc
}).sort_values("porc_importancia", ascending=False, ignore_index=True)
print(f"\nNúmero de features seleccionadas por OMP (k={K_FINAL}): {selected_features.shape[0]}")
print(selected_features.to_string(index=False))
# ---------- Estabilidad de selección (mismos ventana/gap, k=4) ----------
freq = Counter()
for fold, (tr_idx, _val_idx) in enumerate(tscv.split(X_train), start=1):
pipe_f = clone(pipeline).set_params(omp__n_nonzero_coefs=K_FINAL)
pipe_f.fit(X_train.iloc[tr_idx], y_train.iloc[tr_idx])
names_f = strip_remainder_prefix(list(pipe_f.named_steps["prep"].get_feature_names_out()))
coefs_f = pipe_f.named_steps["omp"].coef_
sel_f = list(np.array(names_f)[np.abs(coefs_f) > 0])
freq.update(sel_f)
stab_df = pd.DataFrame({
"feature": list(freq.keys()),
"freq_folds": list(freq.values())
}).assign(porc_folds=lambda d: d["freq_folds"] / tscv.n_splits * 100.0) \
.sort_values(["porc_folds","feature"], ascending=[False, True], ignore_index=True)
print(f"\nEstabilidad de selección (k={K_FINAL}): % de folds en que aparece cada variable")
print(stab_df.to_string(index=False))
# ---------- MODELO FINAL DE PRODUCCIÓN (k=4) ENTRENADO EN TODO EL HISTÓRICO ----------
X_all = pd.concat([X_train, X_test], axis=0)
y_all = pd.concat([y_train, y_test], axis=0)
final_model_k = clone(pipeline).set_params(omp__n_nonzero_coefs=K_FINAL).fit(X_all, y_all)
print(f"\nModelo final 'final_model_k' (k={K_FINAL}) entrenado en TODO el histórico listo para producir predicciones.")
